title: Is this really suited for 2-way ANOVA?

text: My dataset is a biological analysis of a disease which takes 3 months to develop, and we're looking at treatment of that disease. The data is gene expression data and we want to identify changes in gene expression over time. This resembles a 2x2 design in which I have equal sized groups as:                         trtmnt     no-trtmnt     1 month     10           10     3 month     10           10      After thinking about it, I'm not sure if this is really valid for a 2-way ANOVA. The disease develops over the period of 3 months, so the 1 month subjects should look the same regardless of treatment. There are multiple ways to look at this: 1) I could compare trtmnt and no-trtmnt at the 3-month stage. Simple two-group analysis. 2) I could look at 1 month vs 3 month for trtmnt, and then 1 month vs 3 month for no-trtmnt. This would be 2, two-group analyses. But then I would have to compare the two different analyses to look for gene change in common or different. This is where it starts to resemble a 2-way ANOVA, but not really. I appreciate any thoughts. EDIT: I should have pointed out that these are animal studies, so we have 10 animals in each group. Rather than taking samples from each animal twice - once at 1 month and once at 3 months - we have 4 different groups, all of which were used for sampling. I hope that makes sense.


title: How many Folds and Repetitions

text: This is a basic question about training a model. Given the following parameters what would be the criteria I shall use as to to apply a minimum number of folds and repetitions for cross validation of an ensemble of models for regression.               200 rows and 46 columns      300 rows and 53 columns       400 rows and 60 columns        I am unsure where to draw the line as to this. I have read and applied 10 folds and 50 repetitions, 10 folds and 100 repetitions, 10 folds and 5 repetitions and results are mixed. I have read as well that K depends on the percentage of the test sample. It seems that is also data dependant. How can I formulate and apply a minimum number to this two variables.


title: Main effects of 2-way ANOVA - same as 1-way ANOVA on each main effect?

text: I'm working on a 2-way ANOVA analysis with a 2x2 design - two main effects with 2 levels each. In running the analysis, I obviously get the interaction as well as each main effect. Are the results (p and F values) for the main effects the same as if I had simply run a standard ANOVA on each of the main effects independently? That is my interpretation, but when I actually do the test, I get different p and F values.


title: New to the Data Mining business

text: I am new to this data-mining business. I would like for someone to point me in the right direction, plese. In my business - I will have a database with several tables (let's assume MS SQL Server). For the sake of example, I have a 'person' table and an 'purchases' table (each person can make several different purchases). Here is an example of the 'person' table:   * ID   * First Name   * Last Name   * Address   * Date of Birth Here is an example of the 'purchase' table:   * Date of Purchase   * Time of Purchase   * Store Location   * Number of Items Purchased   * Total Purchase Amount What I would like to do is: I would like to perform 'data-mining' in order to find out the "categories", or "groups" of characteristics, that define each 'purchase amount'. For example - I would like to know what is common for all the people who bought between 2,000\$ and 2,500\$ ? Was it that they all came from the same neighbourhood ? Maybe all of them were young & came from the same neighbourhood? Maybe there are 2 distinct different groups that fit the '2000\$-2500\$' spending group ? I am almost sure there are different categories (or distinct fields) that match each price group... one might be old people that come to buy at night, another might be middle-aged people who buy more than 10 items, and another might be middle-aged people who buy one Sundays. From what I understand, this falls under the 'data-mining' region (clustering maybe?) Is this true ? and if so, what are the tools I must use in order to find out these set of parameters that define the different groups for each pricing category? Thank you very much ! John


title: Significant model with no significant betas?

text:
I am running 3 models on 3 subsets of the same data. The set up is as follows:   1. Outcome (DV) is binary categorical   2. Time (IV) is repeated twice (pre and post)   3. Treatement (IV of interest) is binary categorical I am interested to know if at time 2 treatment has had an effect on outcome. I used the lme4 package and used the following R code:               tot.null<-lmer(as.factor(outcome)~Time+(1|id), family=binomial(link='logit'),                  data=df.total)     tot.mod<-lmer(as.factor(outcome)~trt*Time+(Time|id),                   family=binomial(link='logit'), data=df.total)     anova(tot.null,tot.mod)     summary(tot.mod)      **Data head**                  id             trt Time outcome     1   1 peer discussion   -1       1     2   2 peer discussion   -1       1     3   3 peer discussion   -1       0     4   4 peer discussion   -1       1     5   5 peer discussion   -1       1      **str of data**               > str(df.total)     'data.frame':   872 obs. of  4 variables:      $ id     : int  1 2 3 4 5 6 7 8 9 10 ...
          $ trt    : Factor w/ 2 levels "peer discussion",..: 1 1 1 1 1 1 1 1 1 1 ...      $ Time   : num  -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
          $ outcome: num  1 1 1 1 1 1 1 0 1 0 ...      The problem is I get an error messoge on the `tot.mod`:               > tot.mod<-glmer(as.factor(outcome)~trt*Time+(Time|id), family=binomial(link='logit'),                    data=df.total)     Warning message:     In mer_finalize(ans) : false convergence (8)      I think this is the reason the model is significant but none of the predictors are. look at the inflated SEs. **Comparison to the null model and the summary of full model**               > anova(tot.null,tot.mod)     Data: df.total     Models:     tot.null: as.factor(outcome) ~ Time + (1 | id)     tot.mod: as.factor(outcome) ~ trt * Time + (Time | id)              Df    AIC    BIC  logLik  Chisq Chi Df            Pr(>Chisq)         tot.null  3 689.54 703.85 -341.77                                             tot.mod   7 410.67 444.07 -198.34 286.86      4 < 0.00000000000000022 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      > summary(tot.mod)     Generalized linear mixed model fit by the Laplace approximation      Formula: as.factor(outcome) ~ trt2 * Time + (Time | id)         Data: df.total         AIC   BIC logLik deviance      410.7 444.1 -198.3    396.7     Random effects:      Groups Name        Variance Std.Dev. Corr        id     (Intercept)  396.46  19.911                      Time        1441.98  37.973   0.470      Number of obs: 872, groups: id, 436          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)        (Intercept) 10.09866    3.33921   3.024  0.00249 **     trt21        0.01792    5.10796   0.004  0.99720        Time        -0.93753    5.79560  -0.162  0.87149        trt21:Time  -0.84882   10.41073  -0.082  0.93502        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Correlation of Fixed Effects:                (Intr) trt21  Time       trt21      -0.654                   Time        0.558 -0.365            trt21:Time -0.311  0.473 -0.557      What's going on? Why is the model significant but none of the betas? In OLS I know this is an indicator of multi-colinearity among predictors. I don't think that's the reason here. Please help with understanding this problem as well as the error message (I think they may be connected). What are some things I should check for? The other two models from the same data set (`split` on a different grouping variable) had no apparent problems. Thank you in advance. _Using R 2.14.2, lme4 v. 0.999375-42 on a win 7 machine_


title: Mean and Variance of a Set of Sequences

text: I am trying to find the probability distribution (like mean and variance) of a set of sequences which have different length. I came into different measurement models that quantify the distance between difference sequences, but I cannot find something on the probability distribution of sequences. Is there anyway to extract the parameters for the distribution of a set of sequences?


title: Question with R and repeated measures

text: Second time I'm turning to this site, it has helped before. I am looking for some coding help to perform a repeated measures anova in R. I have a total of 4 treatments, each replicated 3 times. They were monitored daily for several months, giving me a value of X each day. So the data is laid out like so:               Day     Trt     Rep     X     1       1       1       0     1       1       2       0     1       1       3       1     1       2       1       0     1       2       2       0     1       2       3       0     1       3       1       0     1       3       2       0     1       3       3       0     1       4       1       1     1       4       2       1     1       4       3       1      What is some example code that would allow me to perform a repeated measures ANOVA on this data?


title: What type of test to use to determine correlation/relationship between two non-continuous varaibles

text: I need to determine if there is any relationship between two count variables. I have 60+ observations for 4 variables and I want to see if any of the pairs of these variables are significantly correlated with one another. Mostly I use R, so forgive me if you're not familiar. I have been using the `cor(...,method="pearson")` and `cor.test()` functions to test each pair, but now I'm not so sure that this is the right approach/test. Would a non-linear regression like `glm(...,family="poisson")` be more appropriate? I started thinking like this because when I looked at a histogram of the counts across my observations, I noticed that there seemed to be a slight tendency for the pink and green variables to go up and down to together. I produced a scatter plot of each of the variables plotted against each of the other variables. I used the tests mention above to try and quantify this relationship and to test weather it was real or just noise. ![enter image description here](http://i.stack.imgur.com/rNFTE.png) ![enter image description here](http://i.stack.imgur.com/iZ3hW.png)


title: Estimate the size of a population being sampled by the number of repeat observations

text: Say I have a population of 50 million unique things, and I take 10 million samples (with replacement)... The first graph is I've attached shows how many times I sample the same "thing", which is relatively rare as the population is larger than my sample. However if my population is only 10 million things, and I take 10 million samples, as the second graph shows I will more often sample the same thing repeated times. My question is - from my frequency table of observations (the data in the bar charts) is it possible to get an estimate of the original population size when it is an unknown? And it would be great if you could provide a pointer to how to go about this in R. ![alt text](http://i.stack.imgur.com/XOEyE.png)


title: Classification after factor analysis

text: I have analysed several dimensions in a survey. Each part of the survey represents a theoretical dimension and is analysed with factorial analysis. I want to use scores from factor analysis to do a classification.   1. The first factors represents a large part of the variance. Can I keep only first factor or do I need to retain all factors?   2. After factor analysis i did a PROMAX rotation which is an oblique rotation. How should i use the output from the PROMAX rotation? If i have taken account of that how do i compute distance with factors correlations matrix ?


title: Post hoc tests in ANCOVA

text: **Question:** What is a good method for conducting post hoc tests of differences between group means after adjusting for the effect of a covariate? **Prototypical example:**   * Four groups, 30 participants per group (e.g., four different clinical psychology populations)   * Dependent Variable is numeric (e.g., intelligence scores)   * Covariate is numeric (e.g., index of socioeconomic status)   * Research questions concern whether any pair of groups are significantly different on the dependent variable after controlling for the covariate **Related Questions** :   * What is the preferred method?   * What implementations are available in R?   * Are there any general references on how a covariate changes procedures for conducting post hoc tests?


title: How does including lags imply losing observations

text: ![enter image description here](http://i.stack.imgur.com/kCFOf.png) My time series professor passed a comment today, that running this model implies losing k+1 observations. I think I have understood how the observations are wasted, but I am notorious for falsely believing that I have understood something(self-diagnosed). A simple explanation would be appreciated!


title: finding probability density function $U=Y-X$

text: suppose joint probability density function of $X,Y$ is $$f(x,y)=\frac{1}{2}\\\ x>0 ,y>0,x+y<2$$how can i find probability density function $U=Y-X$


title: Calendar Visualization

text: I have a two column table consisting of: `date`, `user_id`.   I'd like to visualize a grid with `date` along the `x` axis, and `user_id` along the `y`.   If `date`, `user_id` exists, I'd like to show a grid with a color, say green.   If `date`, `user_id` does not exist, I'd like to show a grid in white.   Bonus points if `date`, `user_id` shows up multiple points resulting in a darker green color. What tool, preferably free and online, is best suited to visualize these data in the manner I've described? Follow up question: which other tools are good at visualizing time series data of user events as a calendar? Something like: https://developers.google.com/chart/interactive/docs/gallery/calendar


title: Thinking about Neural Networks

text: I am confused about the output I am getting for a neural network analysis. What I am trying to do is take in the features from features.csv and use them to create a neural network to predict the variable `Murder.per.100.000`. The first column is a label for each observation so I exclude that. The problem I am having is that `preds` is a vector of `NA`. Some of my features are sparse. I assume I am not dealing with missing values correctly. How should I fix that? When I limit it to non-sparse features, `preds` is made up of all the same value somewhere near 1 which doesn't seem right either.               library('nnet')     dat          = read.table("features.csv", sep=",", header=TRUE)       # remove observation label     dat          = dat[ ,-c(1)]     nData        = dim( dat )[1]     trainInd     = sample( (1:nData), 300 )       # create the training data:     trainData    = dat[ trainInd, ]       # create the test and validation dataset:     notTrainData = dat[ -trainInd, ]     nData        = dim( notTrainData )[1]          valInd   = sample( (1:nData), nData/2 )     testData = notTrainData[-valInd,]     valData  = notTrainData[valInd,]     net      = nnet(Murder.per.100.000 ~ ., data=trainData, size=5, rang=.1,                      decay=5e-4, maxit=200, na.action=na.exclude)     preds    = predict(object=net, newdata=valData)


title: What if probabilities are not equal in the ".632 Rule?"

text: This question is derived from this one about the ".632 Rule." I am writing with particular reference to user603's answer/notation to the extent it simplifies matters. That answer begins with a sample of size $n,$ with replacement, from $n$ distinct items in collection (call) it N. The probability that the $i^{th}$ sample $s_i$ is different from a particular element $m$ of N is then $(1 - 1/n).$ In that answer all elements of N have an equal chance of being randomly drawn. My question is this: suppose instead that in the above question the items to be drawn are such that they are normally distributed. That is, we subdivide the standard normal curve from $Z = -4$ to $Z = 4$ into (say) 100 equal-length subintervals. Each of the 100 items in N has a probability of being drawn that is equal to the area subtended by the curve in its respective interval. My thinking was as follows: The reasoning is similar to that in the linked answer I think. The probability that $s_i \ne m$, with $m$ an element of N, is $P(s_i \neq m) = (1 - F_i)$ in which $F_i$ is the probability of drawing $s_i.$ The probability that a particular element m is in the sample S of size n is $$P(m \in S) = 1 - P(m \notin S) = 1 - \Pi_1^n P(s_i \neq m)$$ $$ = 1 - \Pi_1^n(1 - F_i). $$ A calculation seems to show that as the length of the subintervals gets small the answer converges to the same number as in the first case (probabilities of $s_i$ all equal). This seems counterintuitive (to me) because the construction seems to throw in elements of N which are rare, so I would expect a number smaller than .632. Also, if this is correct, I guess we would have $$ \lim_{n \to \infty} \Pi_1^n(1 - F_i) =\lim (1- 1/n)^n = 1/e, $$ which I don't know to be true or false yet. Edit: If it's true it would probably generalize some. Thanks for any insights.


title: Modeling: Option to cross-validate and predict afterwards

text: This is beyond a R question, But why doesn't it make sense to fit a model, say, a linear discriminant (LDA) model, with _leave-out-one cross validation_ , and afterwards to use this model to predict a new data set? Is this reasoning **right or not** : using such cross validation is a way to measure the performance of the model (instead of testing the model on new test data). If the model performs well, you model it again on the complete training data set and this final model you then use to predict new values.


title: Type I error in t tests and LRT

text: In situations where unequal variance t-test is appropriate, if we use equal variance t-test, type I error would be inflated. I created a likelihood ratio test (LRT) equivalent of unequal variance t-test, but for some reason, it performs poorly (i.e., worse than equal variance t-test). This problem persists even if I increase the sample size to be fairly large (n=100). I don't want it to be just because of the sample size...               NL0 <- function(p,x1,x2){        l1 <- -sum(dnorm(x1,p[1],p[2],log=T))        l2 <- -sum(dnorm(x2,p[1],p[3],log=T))       l1+l2     }          NL1 <- function(p,x1,x2){        l1 <- -sum(dnorm(x1,p[1],p[3],log=T))        l2 <- -sum(dnorm(x2,p[2],p[4],log=T))       l1+l2     }          nsim <- 1000     out1 <- out2 <- out3 <- rep(NA,nsim)     n1 <- n2 <- 20     m1 <- m2 <- 100     s1 <-  1     s2 <- 10          for(i in 1:1000){     x1 <- rnorm(n1,m1,s1)     x2 <- rnorm(n2,m2,s2)          L0 <- optim(c(mean(c(x1,x2)),sd(x1),sd(x2)),NL0,x1=x1,x2=x2)$value         L1 <- optim(c(mean(x1),mean(x2),sd(x1),sd(x2)),NL1,x1=x1,x2=x2)$value          # Test staistic     D <- 2*(L0-L1)     out1[i] <- t.test(x1,x2,var.equal=TRUE)$p.value         out2[i] <- t.test(x1,x2,var.equal=FALSE)$p.value     out3[i] <- 1-pchisq(D,1)     }          mean(out1<0.05)  # equal variance t-test     mean(out2<0.05)  # unequal variance t-test     mean(out3<0.05)  # likelihood ratio test


title: How is this Negative Binomial Random variable used to solve this problem?

text: I was looking at the solution to this problem below and I don't understand how they used a negative binomial R.V. to solve the problem. > A research study is concerned with the side effects of a new drug. The drug > is given to patients, one at a time, until two patients develop side > effects. If the probability of getting a side effect from the drug is $1/6$, > what is the probability that eight patients are needed? The answer that is given to the question is below. > Let $Y$ be the number of patients who do not show side effects. Then $Y$ > follows a negative binomial distribution with $r=2$, $y=6$ and $p=1/6$. > Thus, $$P(Y=6)={(2+6)-1\choose 6}\left(\frac{1}{6}\right)^2 > \left(\frac{5}{6}\right)^6 = 0.0651.$$ Now this is where I'm getting confused. I have my notes on the negative binomial R.V. below: > Suppose that independent trials each having probability $0<p<1$ of being a > success are performed until a total of $r$ successes is accumulated. If we > let $X$ equal the number of trials required then $$ P(X=n)={n-1\choose > r-1}p^r (1-p)^{n-r}.$$   * Why does it seem like $6$ is used for $r$ in one part of the answer but $2$ is used for $r$ in another part of the answer?   * If the value for $Y$ that is used in the answer is supposed to be the same as the value for $n$ in my notes why do they add $2$ to $n$?   * Why is the value for $n-r$ equal to $6$ instead of $4=6-2$? Any clarifying explanations would be really appreciated.


How to prove higher eigenvalues correspond to significant principal components? Is there any explanation or at least intuition for that, if not a proof? I would like to convince myself with some intuition that eigenvalues of covariance matrix correspond to principal components and higher eigenvalues correspond to significant principal components etc.


title: Can a small sample size cause type 1 error?

text: I've learnt that small sample size may lead to insufficient power and type 2 error. However, I have the feeling that small samples just may be generally unreliable and may lead to any kind of result by chance. Is that true?


title: Problem with R code for spectral clustering

text: am writing a simple R script to test the spectral clustering algorithm but for the eigenvalues I don't get them all positive and lambda0 is different from 0. here is my script               Dsqrt<-matrix(0,length(V(g)),length(V(g)))  # The D^-0.5 matrix        for(i in 1:(length(V(g))-1))              Dsqrt[i,i]<-1/sqrt(degree(g,V(g)[i-1]))                     L<-graph.laplacian(g, normalized=TRUE) #Here L[i,i]=1 and L[i,j]= 1/sqrt(d[i]*d[j])        L<-Dsqrt %*% L %*% Dsqrt          Eigns<-eigen(L)      And when I check `Eigns$values` sometimes I get some negative values. Can any one point me my error? Thanks in advance.


title: Does significance test make sense to compare randomised groups at baseline?

text: Many randomised controlled trial (RCT) papers report significance tests on baseline parameters just after/before randomisation to show that the groups are indeed similar. This is often part of a "baseline characteristics" table. However, significance tests measure the probability of getting the observed (or a stronger) difference by chance, aren't they? And if the test is significant we conclude that there is a **true difference because a random difference of that extent would be unlikely**. Does a significance test make any sense after randomisation when we know that any **difference must be due to chance**?


title: Controled Factorial Experiment Design: which name or resource?

text: I'm designing an experiment which has to factors (A,B) with respectively 4 and 5 levels. Since this is quite a large amount, I figured I use a factorial design. However, it is unclear whether the treatment will work at all, at any level or factor. So I designed the following experiment: > All participants perform a test that is a negative control (no > additions/treatment) and then those same participants perform a test with a > randomly assigned level of each factor. This allows me to test if the treatment at any level has any effect (T-test) and then check which combination of factor/level specifically is effective (ANOVA). I've been looking for a proper resource on this, I assume this is a factorial design, but the addition of the control over all participants before the random treatment makes it hard to find the proper title/resource on this design. Can anyone point me in a proper direction?


title: Repeated-measures, crossover analysis using linear mixed model in R

text: I have a data set that I am attempting to analyse in R and I am relatively new to the environment. My full data set contains 7 subjects (represented by Subject), that all receive 3 treatments (environmental conditions, represented by Altitude) and are measured 10 times within each treatment. Each participant received each treatment in a different order, represented The measurements are power (from a treadmill, represented by Power) during repeated-sprint efforts (represented by Sprint). An example: ![enter image description here](http://i.stack.imgur.com/lHITL.jpg) I also have characteristics about each subject that may influence their power effort, including weight, capacity tests and age. I am interested in the effect of each Environment on the Sprint results, plus the effect of order on the Power. I believe a linear mixed model, with subject as a random effect, is an appropriate tool to investigate my dataset. I have attempted this using the following line: alt.model = lmer(Power ~ Sprint + Altitude + VO2max + (1|Subject), data=ALTMM) However, I don't think this accounts for each individual sprint. How do I represent this? Thank you.


title: What approach to use for attrition analysis?

text: I am trying to perform an attrition analysis on a company with an average size fo around 250 to 300. I have monthly attrition data for the last 30 months or so. Now if i want to go for a predictive approach where i can predict the attrition for next few months, say 3 months, which one would be the best approach. After some research I am juggling between poisson's distribution and time series analysis. Which approach is the better way to go ?? and why ?


title: Fast integration of a posterior distribution

text: I wish to infer the posterior distribution on the probability of success $\theta$ in some binomial process, the twist being that I know that $\theta$ lies in the interval [0.5, 1]. The trouble is that a Beta distribution that is supported on [0.5, 1] is no longer conjugate. In particular, the prior looks like this ($p$ is the pdf): $p(\theta) \propto {(\theta-0.5)}^{\alpha-1} {(1-\theta)}^{\beta-1}$ The likelihood, given $x$ successes in $n$ trials, is the usual binomial likelihood: $p(n,x | \theta) \propto \theta^x {(1-\theta)}^{n-x}$ The posterior, as it stands, looks like this: $p( \theta | n, x) \propto \theta^x {(\theta-0.5)}^{\alpha-1} {(1-\theta)}^{n-x+\beta-1}$ I would like to be able to make CDF calculations (i.e. compute integrals) of the posterior efficiently. (numerical integration, e.g., may be too slow for my purposes). I am unable to figure out how to do this. Much appreciate any help the forum can provide.


title: What subject is more useful for statistics students

text: I am deciding which of the following two subjects to take as one elective in the coming term.   1. Complex Analysis   2. Theory of Partial Differential Equations Ideally I'd like to take both. However, I will have other things to do. Therefore, I want to take something that is direct relevant to statistics and probability. Any suggestion? Thank you!


title: To get the overall significant p value, am I on the right track for LMM?

text: According to this document: http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf I can get the overall p-value of fix effects by comparing models that I would like to know to null model. In my case I would like to know the effect of three-way interaction. So first, I run the null model:               lmer86 <- lmer(IntensityMax ~ (1|item) + (1+vowel3|speaker) + sex + vowel3 + Language, data=data1.frame, REML=FALSE, na.action=na.omit)      Then I run the model that I would like to test:               lmer81 <- lmer(IntensityMax ~ (1|item) + (1+vowel3|speaker) + sex*vowel3*Language, data=data1.frame, REML=FALSE, na.action=na.omit)      And finally, I compare these models using ANOVA function:               anova(lmer81,lmer86)      And this is result:               Data: data1.frame     Models:     lmer86: IntensityMax ~ (1 | item) + (1 + vowel3 | speaker) + sex + vowel3 +      lmer86:     Language     lmer81: IntensityMax ~ (1 | item) + (1 + vowel3 | speaker) + sex * vowel3 *      lmer81:     Language            Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)         lmer86 14 13349 13432 -6660.5    13321                                  lmer81 26 13319 13474 -6633.6    13267 53.813     12  2.951e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Then I know that three-way interaction is an important effect. Am I on the right track?


title: Calculating hypergeometric probability with approximations

text: I'm essentially trying to reproduce the `dhyper` function in R using   `!n ~ gamma(n+1)`. Why does this approach not work for large values of the population?               # (n choose k)     fastChoose <- function(parm_n, parm_k) {        p1 = lgamma(parm_n+1)        p2 = lgamma(parm_n - parm_k +1)        p3 = lgamma(parm_k + 1)        result = p1 - p2 - p3        return (exp(result))     }          # P(X=k) = (K choose k) * (N-K choose K-k) / (N choose n)     myDHyper <- function(parm_k, parm_K, parm_N, parm_n){        return (fastChoose(parm_K,parm_k) * fastChoose(parm_N-parm_K,parm_K - parm_k) /     fastChoose(parm_N,parm_n))      }          BIGN<-15000     whiteBallsInPop <- 50     iSampled<- 200     allPointMass<-c()     allPointMassMine <- c()     for(i in seq(0,whiteBallsInPop)){            fromDHyper <- dhyper(i, whiteBallsInPop, BIGN-whiteBallsInPop, iSampled, log = FALSE)       allPointMass <- c(allPointMass, fromDHyper)            fromMyDHyper <- myDHyper(i,whiteBallsInPop,BIGN,iSampled)       allPointMassMine <- c(allPointMassMine,fromMyDHyper )          }           sum(allPointMass)     sum(allPointMassMine)      Keeping everything on log scale and taking log(p) doesn't seem to work either:               #log( P(X=x))     myDHyperOverflowFixed <- function(k, K, N, n){       # P(X=k) = (K choose k) * (N-K choose K-k) / (N choose n)       #        = K!/k!(K-k)! * (N-K)!/(K-k)!(N-K-K+k)! /  N!/n!(N-n)!             log_p <- lgamma(K + 1) +             lgamma(N - K + 1) +             lgamma(n + 1) +             lgamma(N - n + 1) -              lgamma(k + 1) -               lgamma(K - k + 1) -               lgamma(K - k + 1) -               lgamma(N - K - K + k + 1) -              lgamma(N + 1)           return (log_p)     }


title: Confusion related to online bounds for bayesian algorithms

text: I was reading this paper Online Bounds for Bayesian algorithms and they had some derivations. I didn't get how they arrived to the conclusion ![enter image description here](http://i.stack.imgur.com/tIGtF.png) ![enter image description here](http://i.stack.imgur.com/YjE4C.png) I didn't get how equation 3 was derived any suggestions guys?


title: Logistic vs. Linear regression for recovery rate modeling

text: i am trying to model recovery rates in my data that are in the range of (-.1,+1). Around 8% of my observations are negative too. Broadly, below is the dist.: RR<0: 8% RR=0: 30% RR>0 and RR<0.95: 53% RR>0.95: 8% what would be the best way to model recovery rate: linear or logistic events and trial?


title: strucchange package on ARIMA model

text: Is there a way to use strucchange package in R on ARIMA models? I haven't been able to find any. Thanks a lot.


title: What's the difference between principal components analysis and multidimensional scaling?

text: How are PCA and classical MDS different? How about MDS versus nonmetric-MDS? Is there a time when you would prefer one over the other? How do the interpretations differ?


title: Mixed variable, joint distribution, How do we know which one is continuous distribution, which one is discrete

text: If we have one continuous r.v. $x$ and a discrete r.v. $y$ which takes one of the two values $y_1$ and $y_2$. Let's say we know the prior probabilities $P(y_1)$ and $P(y_2)$. From Bayes theorem we know $p(y_i,x) = P(y_i|x)\cdot p(x)$ $P$ corresponds to prob. mass function, $p$ corresponds to prob. density function. On LHS we have a probability value coming from a continuous distribution. On RHS we have product of two probabilities one coming from a discrete distribution, another coming from a continuous distribution. My question is - Why is first term on RHS $P(y_i|x)$ and not $p(y_i|x)$ (Since it's possible that, given $x$, the probability that $y$ will take a value of $y_i$, might smoothly vary from $0$ to $1$)? Are there different ways of writing Bayes theorem that will have $p(y_i|x)$ instead of $P(y_i|x)$?


title: The efficiency of Decision Tree

text: When we use top down approach to induce a decision tree, we need to use some kind of splitting criteria to choose the _splitting feature_ and _splitting value_ at a certain internal node. The criteria could be:   * IG   * IGR   * Gini impurity But the feature selection seems **_very computationally intensive_**. Because we need to enumerate all the existing features and go through all its possible splitting values. And for each of the feature/splitting value, we need to calculate the criteria value and memorize it, and select the best value. Am I correct on this? Is there anything we can do to improve it?


title: After install a package in R do we need to require all the time?

text: I have installed rjson in R. When I try to call a function from this package like toJSON it give a error like undefined function. But if I use require("rJSON") before to call the function it works well. SO do I need to do that all the time? Thanks


title: How to test whether ICCs from different groups are significantly different?

text: I need to compare two generalizability (G) coefficients for data that are from two separate populations. G coefficients are a type of intraclass correlation coefficient (ICC). The literature on statistically comparing G coefficients is sparse by all accounts, so I thought that the best I might do was find out the recommendation for statistically comparing ICCs. I've considered Fisher's Z, but I know ICCs are not normally distributed.


title: multiple choice data simple logistic regression or multinomial logistic regression

text: I have a survey question where the respondent can check one choice or two choice maximum. Question looks like this: > What is the more important characteristic when you buy chocolate? > >   * Sweetness (characteristic 1) >   * Packaging (characteristic 2) >   * Price (characteristic 3) >   * ... (characteristic 4) >   * ... (characteristic 5) >   * ... (characteristic 6) > In my data, I have 6 column (1 for each characteristic) with the value 1 if the respondent checked the ad-hoc column and 0 otherwise. I would like to perform a regression where the outcome is the question and the predictors the socio demographic criteria (some of them are numerical such as the age and other categorical such as the gender). Now I wonder how I should do. Should I do a simple logistic regression 6 times with same predictors but each time different outcome? >   * glm(caracteristic1 ~ age+gender) >   * glm(caracteristic2 ~ age+gender) >   * glm(caracteristic3 ~ age+gender) >   * glm(caracteristic4 ~ age+gender) >   * glm(caracteristic5 ~ age+gender) >   * glm(caracteristic6 ~ age+gender) > Or is it possible to use another method such as the multinomial logistic regression. If yes which R packages do you recommend me to use?


title: (R) All data going into 1st row of new column of data frame

text: I've got a data frame, 1447 obs of 165 variables, and am looking to add columns by processing existing ones. I've got:               resvar <- 111       # column number of response variable (e.g. CPUE) in "samples" (the data frame)      and am looking to log1p that column, 111, row by row, to a new column, called "grv" at the end of samples:               samples$grv <- log1p(samples[resvar])      But it puts the output into column grv, row 1, as ``c(item1,item2,...item1447)`               class(log1p(samples[resvar])) = data.frame      Does anyone know what's going on & why it's not working? If I do               dummy <- log1p(samples[resvar])      `dummy` is a data.frame with row name the same as the row name from resvar, column 111's row name. Is this the problem?


title: Comparing the mean of two non-random sub-groups

text: Students are assigned to group A or B based on two aspects: a grade from an exam and an evaluation of their teacher. Students with both a high grade and a positive evaluation get into group A, other students in group B. I know that the grade is much more influential than the teacher’s evaluation. I have the exam grades of 1,000 students. I also know whether students are assigned to group A or B. I want to know whether the average grade of students in group A is higher than the average grade in group B. I am thinking about this regression (which I believe is equivalent to a direct t-test): $\text{Grade}_i=\alpha+\beta*\text{DumA}_i+\epsilon_i$ where $\text{Grade}_i$ is the grade of student i. $\text{DumA}_i$ equals 1 if student i is in group A and 0 otherwise. Do students in group A have a higher grade if Alfa is significant according to a standard t-test? Note that students in group A tend to have a higher grade since assignment is mostly based on the grade. It seems to me that the dummy is endogenous and thus the estimate of Alfa is incorrect. Is my interpretation correct? Does this mean that it is not possible to compare the mean of the two groups? Are there any alternatives? Any thoughts are very much appreciated!


title: Statistical hypothesis test on euclidean spaces

text: On Lehamman and Romano's Testing Statistical Hypothesis, 2nd edition, there is a problem that says "If the sample space is euclidean and $P_0$ and $P_1$ have densities with respect to the lebesgue measure, there exists a nonrandomized most powerful test at every significance level $\alpha$." Why is it not a direct consequence of the Neyman-Pearson theorem?? The hypothesis are both simple and the $P_0$ and $P_1$ have densities!


title: Kaczmarz's projection algorithm

text: I am trying to understand the derivation of the following formulas given in my lecture (sadly without any further explanation). It says that each new prediction error is of the form $e(r+1)=y(r+1)-h^T(r+1)\cdot \pi(r)$. The problem yield the following minimization problem $J(\pi)=\frac{1}{2}\cdot (\pi(r+1)-\pi(r))^T\cdot (\pi(r+1)-\pi(r)) + \mu \cdot[y(r+1)-h^T(r+1)\cdot \pi(r+1)]$ In the lecture it was said that the necessary conditions for the min. solution are: $\dfrac{\partial J}{\partial \pi(r+1)} = \pi(r+1)-\pi(r)-\mu \cdot h^T(r+1) = 0$ $\dfrac{\partial J}{\partial \mu} = y(r+1) - h^T(r+1)\cdot \pi (r+1) = 0$ So far so good. What I don't understand is how to get the following obtained solution "by solving the set of linear equations": $\pi(r+1) = \pi(r) + \dfrac{h(r+1)}{h^T(r+1)\cdot h(r+1)}\cdot [y(r+1)-h^T(r+1)\cdot \pi(r)]$ How do I get this solution? Any help is highly appreciated


title: Full covariance matrix of multivariate mixture Gaussian distribution

text: I want to estimate the full covariance matrix of a multivariate-Mixture Gausian distributed vector. I have a multivariate-Mixture Gausian distribution for a vector $\mathbf{r}=[r_1,..,r_T]^T$ on the form \begin{align} p(\mathbf{r})&=\sum_{\boldsymbol{\pi} \in \Omega_{\boldsymbol{\pi}}^T} p(\mathbf{r|\boldsymbol \pi})\prod_{t=1}^T p(\pi_t),\\\ \end{align} The latent categorical variable $\boldsymbol{\pi}=[\pi_1,...,\pi_T]$, takes values from a discrete state space such that $\pi_t \in \Omega_{ \pi} : \\{1,...,L\\}$ and therefore $\boldsymbol{\pi} \in \Omega_{\pi}^T$. where I have decoupled the variance matrix in the following manner for simplicity: \begin{align} p(\mathbf{r}|\boldsymbol{\pi})=&N_T\left( \begin{bmatrix} \mu_{r_1|\pi_1}\\\ .\\\ .\\\ .\\\ \mu_{r_T|\pi_T} \end{bmatrix}, \begin{bmatrix} \sigma_{r_1|\pi_1}^2 & & & &0\\\ & .& & & \\\ & &. & & \\\ & & &.&\\\0& & & & \sigma_{r_T|\pi_T}^2 \end{bmatrix} \begin{bmatrix} 1 & \vartheta(1) & .& .& &\vartheta(T-1)\\\ \vartheta(1) & 1& & & &\\\ .& &. & & & \\\ .& & &.& & \\\\\vartheta(T-1) & & & & & 1 \end{bmatrix} \right)\nonumber \\\ \nonumber \\\ =&N_T\left(\boldsymbol{\mu_{r|\pi}},\boldsymbol{\Sigma_{r|\pi}^{\sigma}} \boldsymbol{\Sigma_{r|\pi}^{\vartheta} }\right)\\\ =&N_T\left(\boldsymbol{\mu_{r|\pi}},\boldsymbol{\Sigma_{r|\pi}}\right)\\\ \end{align} How do I find the full covariance matrix of $\mathbf{r}?$ I have tried \begin{align} var(\mathbf{r}) & = E ( var(\mathbf{r}\mid \boldsymbol{\pi})) + var (E(\mathbf{r} \mid \boldsymbol{\pi})) \\\\[12pt] & = \sum_{\boldsymbol{\pi} \in \Omega^T} p(\boldsymbol{\pi}) \boldsymbol{\Sigma_{r|\pi}} + \sum_{\boldsymbol{\pi} \in \Omega^T} p(\boldsymbol{\pi}) (\mu_{r|\pi}-\bar{\mu_r})(\mu_{r|\pi}-\bar{\mu_r})^T \end{align} But here you have to sum over L^T terms. Any tips on how to estimate the full covariance matrix of $\mathbf{r}$ is very welcomed.


title: Which one consider when dealing with dummy variables: Exact Sig. (2-sided) or Exact Sig. (1-sided)?

text: I'm running a cross tab of one categorical variable against a dummy variable. This dummy variable indicates the presence or absence of a company characteristic(if the company is industrial it should be 1 and 0 otherwise). I this case, when running the exact test, shoul I consider the Exact Sig. 1 sided or the 2 sided? Thanks


title: Estimating Probabilities from Stratified Table

text: Say I have data stratified by smoking status, where $Y=1$ is has lung cancer, $Y=0$ doesn't have lung cancer, and finally $D=1$ drinks alcohol, and $D=0$ doesn't drink alcohol. How do I estimate $P(Y=1| D=1, S=1)$? It follows this formula I believe(http://math.stackexchange.com/questions/419408/how-to-compute- conditonal-probability-on-multiple-conditions-on-the-right-side) First I look at Smokers chart- and under $Y=1$ & $D=1$ I see the number is $24$. Is this over the probability that $D=1$ and $S=1$, so now it is $24/(24+776)=24/800$? And finally $P(Y=1|D=0, S=0)= 21/(21+2079)=21/2100$? **Smokers ($S=1$):**                   |  Y=1    Y=0      ------------------         D=1 |   24    776     D=0 |    6    194      **Non-Smokers ($S=0$):**                   |  Y=1    Y=0     -----------------     D=1 |    9    891     D=0 |   21   2079


title: Cantelli's inequality proof

text: I am trying to prove the following inequality: EDIT: Almost immediately after I posted this question, I discovered that the inequality I am being asked to prove is called Cantelli's inequality. When I wrote this up, I didn't realize this particular inequality had a name. I have found multiple proofs through Google, so I don't strictly speaking need a solution anymore. However, I am keeping this question up because none of the proofs I have found involve invoking the fact that $t=E(t-X)\leq E[(t-X)\mathbb{I}_{X<t}]$, as was originally intended. For $t \geq 0$, $\mathbb{P}(X-E(X)\geq t)\leq \frac{V(X)}{V(X)+t^2}$ Our professor gave us the following "hints" for working this out: "First work out the problem assuming $E(X)=0$ then use the fact that $t=E(t-X)\leq E[(t-X)\mathbb{I}_{X<t}]$." EDIT: To be clear, in my notation, $\mathbb{I}$ refers to the indicator function. The first part is pretty simple. It is basically a variation of the proof for Markov's or Chebychev's inequality. I did it out as follows: $V(X)=\int_{-\infty}^{\infty}(x-E(X))^2f(x)dx$ (I know that, properly speaking, we should replace $x$ with, say, $u$ and $f(x)$ with $f_x(u)$ when evaluating an integral. To be honest, though, I find that notation/convention to be unnecessarily confusing and not terribly transparent, so I am sticking with my more informal notation.) If we assume $E(X)=0$, then the above simplifies to $V(X)=\int_{-\infty}^{\infty}x^2f(x)dx$ For brevity's sake, I will skip some steps, but it is easy to show then that $V(X)\geq t^2 P(X>t)$, or rather $P(X>t)\leq \frac{V(X)}{t^2}$. Since $E(X)=0$, we can replace the $X$ on the left-hand side of the latter with $X-E(X)$. Here is where I am having trouble moving forward. I don't understand how to go about using the fact that $t=E(t-X)\leq E[(t-x)\mathbb{I}_{X<t}]$. Again, since $E(X)=0$, we can substitute in $t-E(X)$ for $t$. This is equivalent to $E(t-X)$. Then, we can rewrite the $t^2$ in the denominator at the right-hand side of the inequality as $[E(t-X)]^2$, which since the middle term drops out simplifies to $t^2-[E(X)]^2$. But I don't see where I can go from here, either. Though you can further rewrite this as $t^2+V(X)-E(X^2)$, which at least gets me the $V(X)+t^2$ term in the right place. Clearly I am missing something, here, related to $E(t-X) \leq E[(t-X)\mathbb{I}_{X<t}]$, but I quite frankly just have no idea what to do with this term. I understand conceptually what this term is telling me. Intuitively, the expected value of $t-X$ is going to be smaller than the same quantity if $X$ is restricted to being strictly less than $t$; that is, the former term is likely to be negative, while the latter must be positive. But I don't see how I can use this fact in the proof. I tried "distributing" on the inside to simplify ... $E[(t-X)\mathbb{I}_{X<t}]=E[t\mathbb{I}_{X<t} - X\mathbb{I}_{X<t}]=tP(X<t)-?$ But am not sure how to evaluate $E(X\mathbb{I}_{X<t}]$. Anyone have an idea or a hint?


title: Clustering data for occurrence

text: I have a set of data representing nodes and how often they have been involved with each other. I've processed this into a table containing the nodes on X and Y with the data being the number of interactions e.g.                  A    B    C    D     A  2    2    0    0     B  2    2    0    0     C  0    0    2    1     D  0    0    1    2      So the data can be looked at either way. I would now like to perform some clustering of these to find clusters of related nodes through their level of interaction, as well as some visualisation. I have R which is ideal but am happy to use something else if that's easier. All the for-noddy clustering I see doesn't have this type of data: it would say have a set of nodes and for each a set of variables and clustering is for similarity in those variables - a different problem than I'm attempting.


title: What makes the canonical link function special in GLMs?

text: Why is the canonical link function used so frequently with GLMs? What makes it "natural"? Is there any reason to think that, $Q(\theta _i)$ (where $Q$ is the canonical link function, and $\theta _i$ is the parameter of interest) is better described by a linear combination of predictor variables than some other function of $\theta_i$. That is, is there any reason to believe that: $$Q(\theta _i) = \sum_j \beta _j x_{ij}$$ is superior to: $$f(\theta _i) = \sum_j \beta _j x_{ij}$$ where $Q \neq f$ If not, is there anything else that makes $Q$ better than $f$? The text book I am using just mentions what a canonical link function is, and makes use of it (pretty much exclusively), but does not explain what distinguishes it from any arbitrary link function.


title: Using a binomial approach to calculate significance of multiple significance tests

text: I have 10 data sets with normal distributions and carrying out significance tests on each one. With this many data sets, there is a good chance that some will be significant by chance. I am going to try the Bonferroni method to solve this problem but another approach that occurred to me is a binomial one. i.e. if the p value for each individual significance is 0.01 then the probability of getting a significant result follows B(10, 0.01). You can then generate the probability of any number of significant results. Is this method "legit"? I can't see anything wrong with it and it has been mentioned in another thread but does not seem to be very common and woud like to know if there is an issue with using it.


title: How to interpret if my sample statistic is way out in the tail of the bootstrap distribution

text: I use bootstrapping to generate the distribution / histogram of my sample statistic and find out that the value of my real sample statistic is way out in the tail. What does this mean? Does it mean that my sample is unrepresentative of the population? But if it is indeed unrepresentative of the population, then wouldn't bootstrapping give the incorrect sampling distribution in the first place? This is a "what if" question of mine, so if it is somehow impossible for my real sample statistic to be way out in the tail, I'd appreciate an explanation about why as well.


title: How is the energy function of a Gaussian RBM derived?

text: I think it must have been come through the PDF of a Gaussian distribution. Could somebody elaborate?


title: What does drawing sample using Metropolis-Hastings algorithm mean?

text: I am confused with the word " **draw samples from any probability distribution P(x)** ", mean I apologize for my ignorance, but, drawing sample as i understand, is for example, tossing a coin and writing whether the outcome was a head or a tail, x1 {1,0}, whether a subject was obese, normal or underweight , x2 {1=Obese, 2=Normal, 3=Underweight }, am I wrong ? So please help me understand what " **draw samples from any probability distribution P(x)** ", mean ?


title: How to test difference between weighted and non-weighted slope in R?

text: I have a dataframe with the following column: Year, Temperature, Num.Species which are all numeric variables. Each row represent the number of species observed a certain year and the mean temperature that year. I want to test if the regression slope of Year~Temperature, differs from the regression slope of Year~Temperature when using Num.Species as a weighting factor. From what I understood so far, it is possible to do this with an ANCOVA. But there is no categorical variable in my data. Any help would be appreciated!


title: 3 way mixed ANOVA & pairwise comparisons effect size

text: im in desperate need of some help for my data analysis assignment! I have a 3 way mixed anova- where I have 2 sig main effects, one significant two way interaction and the sig three way interaction Ive been told to not describe the lower order effects if the three way is significant so I haven't ( but have put inferentials in table) After getting the sig 3 way I ran a post hoc test to compare the means, using the EMMEANS command ( as per the handout I got). As its anova it gives partial n2, I know how to change the main effects and interaction into eta n2 (using the within effects table and between effects table obviously depending on what sort of factor it is) But with the pairwise comparisons which I got from doing the EMMEANS command ' EMMEANS = TABLES (type*wggroup*match) COMPARE (match) /EMMEANS = TABLES (wggroup*type) COMPARE (type) how do I change the partial n2 in the multivariate table to eta n2?? as I now have the inferentials in my text for the pairwise comparisons but am unable to convert the effect size to eta n2 from the multivariate table to interpret it correctly the only tables I am given from EMMEANS is the pairwise comparisons ( which doesn't have the inferentials) and then below that is the multivariate table ( with inferentials) thanks!


title: Normalize sample data for clustering

text: I have three types of summary score, $a, b$ and $c$ for 200 samples. $S1, S2, S3,..., S200$ $a_{s1}, a_{s2}, ..., a_{s200}$ $b_{s1}, b_{s2}, ..., b_{s200}$ $c_{s1}, c_{s2}, ..., c_{s200}$ Each of these scores is essentially the number of times that some event occurs in the data of each sample. I wish to find groups of these samples by correlation-based clustering. However, the scales for each of these scores are very different: Summary of $a$:               Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      2.0   36.0   55.0   52.5   69.0  139.0       Summary of $b$:               Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      8.0   99.5   285.0   292.7   737.5  2624.0       Summary of $c$:               Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      40.0    111.0   176.0   300.4   554.5   779.0       Should I have to normalize the scores? If so, how? NB. I want to make no assumptions about the dependence or independence between these types of events and hence between these summary scores. UPDATE: So, I've decided to try clustering with Euclidean. In order to get normalized and transformed data, I'm doing the following: 1\. test whether scores fit a normal distribution with Shapiro test   1. if not, then do a boxcox transformation using $\lambda$ from a boxcoxfit   2. get z-score for each score   3. cluster with euclidean distance measure Does this seem reasonable?


title: Relationship between dimentionality reduction and clustering algorithms

text: I've got bit confused about dimensionality reduction and clustering . whether all clustering algorithms (k-means, affinity propagation, spectral clustering,...) do kind of dimensionality reduction ?


title: Proper use of the coefficient of variation

text: I have a variety of samples, each with a different standard deviation and mean. The coefficient of variation $CV$ = ${\sigma} / {\mu}$ defines the amount of variation in a population or sample around its mean. Is it meaningful to then use $1/CV$ or some other variation on this as a weight so that, using the example above, $\mu / CV$ is the new mean for each sample? Is this mean at all more accurate for comparing the samples? By comparing the samples I would like to determine (eyeball) whether, adjusted for standard deviation, the means of the samples differ significantly. I'm almost sure this is an abuse of the coefficient of variation. Are there other less silly methods I should be looking at to adjust for differences in standard deviation among samples?


title: Prediction/generalisation error derivation

text: In section 7.3 of "The Elements of Statistical Learning" (Hastie, Tibshirani, Friedman 2009) one can find the following: ![http://i.stack.imgur.com/IKbnD.jpg](http://i.stack.imgur.com/LIMm4.jpg) **I am at a loss at how to get from the first to the second line, and how to get from the second to the third line in the derivation:**   1. In line 2, I get the first and second term, but the third appears out of nowhere: I only get $E\left[\hat{f}^2\right]$ for the third term when I expand the first line.   2. In line 3, how does the variance follow from line 2? Is this a typo in line 2? In my book, $\operatorname{Var}\left(\hat{f}\right) = E\left[\left(\hat{f}-E\hat{f}\right)^2\right]$, not $\operatorname{Var}\left(\hat{f}\right)=\left(E\left[\hat{f}-E\hat{f}\right]\right)^2$ as suggested.  Can anybody enlighten me?


title: Correlation not significant - interpretation?

text: Sorry that the question may sound stupid, but I'm really confused right now. For getting to know statistics better, I replicated an existing study (same questions) and did the survey (n=66). I calculated correlations (Pearson bivariate, 2 sided) and found statistically significant results. Unluckily, some of my results are insignificant where the first study (which I replicated) found significant results. What can I say about that (besides that I have to reject the hypothesis and I would need a larger sample size to get significant results)? Is the first study wrong?


title: Decaying effect of IV over DV. How to analyze?

text: I need to prove my hypothesis: The relative impact of the direct ties on the project outcome decreases as the direct ties network grows. So I have IV (network size) and DV (outcome). I have proved by regression that IV is a significant predictor over the DV. However, how can I prove that the effect decays as the IV grows?


title: SE of RMSE in R

text: I have crossvalidated my models and measured RMSE between the modelled values and reality:               RMSE <- function(err) sqrt(mean(err^2))     RMSE(predicted - reality)      I am going to choose model with lowest RMSE. But I would like to know what difference in RMSE is still significant. I guess I could do some F-test comparisons, but as I have tens of models I don't like the idea of pairwise comparisons. This leads me to question ### How can my RMSE function be extended to compute SE of RMSE? How to do this in R? I guess this will be somehow based on Chi-square distribution parametrized by `length(err)` and maybe somehow scaled by `sd(err)` (?), but I don't know how to make this statistically correctly in R.


title: How to get observations from residuals in an ARIMA model?

text: If we have residuals of an ARIMA(p,d,q) with known parameters, how can we retrieve the original observations of the time series?


title: Does averaging across multiple time series render higher Nyquist frequency?

text: Suppose I have three time series with the same length and the same temporal interval. If I upsample each of the three time series through interpolation, that will cause aliasing. However, if I'm interested in the information embedded in the average of the three time series, is such an upsampling justified? In other words, does the above upsamping and then averaging among three time series increase the Nyquist frequency by three times? Thanks!


title: What are the differences between Factor Analysis and Principal Component Analysis

text: It seems that a number of the statistical packages that I use wrap these two concepts together. However, I'm wondering if there are different assumptions or data 'formalities' that must be true to use one over the other. A real example would be incredibly useful.


title: Factor analysis using principal component method

text: Is factor analysis using principal component method the same as principal component analysis?


title: PCA vs. Factor Analysis

text: When should we use PCA over factor analysis? Aren't they essentially the same thing except that factor analysis is modeling observed variables as linear combinations of unobserved factors? Whereas PCA is modeling components as linear combinations of observed variables?


title: What statistical test is best for my biological problem?

text: I have a dataset that has abundance of different bacteria in different mice type. Let’s say mice are called: sw/c57, sw/sw, c57/c57 and c57/sw. The fecal samples of these mice have been analyzed and the number of present bacteria in their fecal samples is available. Now I want to find those bacteria that are significantly different in any of the mice type. E.g. Prevotela is significantly different in c57/c57 comparing to other mice. An example of my dataset:                          sw/sw    sw/c57    c57/c57    c57/sw     Roseburia  500      141       0          589     Prevotela  300      131       237        288      Each mice type has 10 replicates. What kind of statistical test I should use to analyze my data? Thanks!


title: R package ltm: How to manipulate title on item response category characteristic curve plot

text: I have conducted an IRT analysis using package ltm in R 2.10.               fit1 <- grm(data)      When producing plots of item response characteristic curves (ICCs),               plot(fit1)      the text "Item response characteristic curves - Item:NAME" is printed as the main title of each plot. NAME is substituted with the actual name of the items in the dataset. A common practice is to present all the individual ICC plots on a single page, and then all that is needed in the title of each plot is the NAME (or at the most "Item: NAME"). I have not been able to achieve this, if I create a vector with the actual item names (e.g. v1, v2, v3...Vn) they are all used as the title in each of the plots. What I would like to achieve is to have the title of each plot simply read "Item: NAME" where NAME is replaced/pasted with the names of the actual variables in my data. The variable names themselves are available in the coefficients `fit1` object, and can be retrieved with `names(fit1$coefficients)`. Anyone have a suggestion for how to achieve this? Reproducible sample (thanks @Chi)               dat <- Science[c(1,3,4,7)]     fit1 <- grm(dat)     par(mfrow = c(2,2))     plot(fit1)      This plots each of the four ICCs on the same device, but all titles contain the abovementioned text that I wish to remove. It also contains the text "Item: NAME" (i.e. comfort, environment, work, benefit) which I would like to retain.


title: correlation of two sums of random variables

text: Imagine two random variables $X$ and $Y$ which are correlated with $\rho = 1$. Both have a mean of $100$ and a standard deviation of $40$. Two other random variables $U$ and $V$ are correlated at $\rho=0.8$. Both have a mean of $0$ and standard deviation of $20$. Now, I wonder if there is a formula to compute the correlation of $A = X+U$ and $B=Y+V$? $\text{cor}(X,U)=\text{cor}(X,V)=\text{cor}(Y,U)=\text{cor}(Y,V)=0$ Any ideas on how to compute $\text{cor}(A,B)$? Thanks in advance


title: Prove that $E[g(X)] = \int_{-\infty}^{\infty}G(t)\phi(t) dt$

text: Let $X$ denote a real-valued random variable with characteristic function $\phi$. Suppose that $g$ is a real-valued function on $\mathbb{R}$ that has the representation $\hspace{25mm}g(x) = \int_{-\infty}^{\infty}G(t)\exp(itx)\,dt$ for some $G$ satisfying $\hspace{25mm}\int_{-\infty}^{\infty}|G(t)|dt < \infty$ Prove that $\hspace{15mm}$$E[g(X)] = \int_{-\infty}^{\infty}G(t)\phi(t) dt$ Attempt: $\hspace{15mm}$ $E[g(X)] = \int_{-\infty}^{\infty}G(t)\exp(itX)\cdot p(x) dx$ Since $\int_{-\infty}^{\infty}\exp(itX)\cdot p(x) dx$ = $E[\exp(itX)] = \phi(t)$, we then have $\hspace{15mm}$$E[g(X)] = \int_{-\infty}^{\infty}G(t)\phi(t) dt$


title: Can I use Synthetic Control Method for Comparative Case Studies with survey data?

text: I'd like to assess the impact of an upcoming policy implementation, as measured by changes in questionnaire response to a Likert-scale question. I understand I could use a difference-in-difference approach. However, in my situation there is no single obvious comparison, non-treated population. I think I'd like to use the "Synthetic Control Method for Comparative Case Studies" as described by Abadie et al and implemented as Synth in R.   * Alberto Abadie, Alexis Diamond, Jens Hainmueller. Journal of the American Statistical Association. June 1, 2010, 105(490): 493-505. doi:10.1198/jasa.2009.ap08746. full text As summarized in the R help for synth: > synth estimates the effect of an intervention of interest by comparing the > evolution of an aggregate outcome for a unit affected by the intervention to > the evolution of the same aggregate outcome for a synthetic control group. > > synth constructs this synthetic control group by searching for a weighted > combination of control units chosen to approximate the unit affected by the > intervention in terms of the outcome predictors. The evolution of the > outcome for the resulting synthetic control group is an estimate of the > counterfactual of what would have been observed for the affected unit in the > absence of the intervention. [..] the synth function routinely searches for > the set of weights that generate the best fitting convex combination of the > control units. In other words, the predictor weight matrix V is chosen among > all positive definite diagonal matrices such that MSPE is minimized for the > pre-intervention period. See also the useful summary by Srikant Vadali in answers below. Is this method appropriate for survey/sampled data? Is there anything I need to do differently, or just use my Likert-response mean as the dependent variable? Any suggestions about how I'd power such a beast? Thank you!


title: How to estimate similarity between several probability distributions?

text: I have several set of probability distributions. I want to reliably estimate consistency across distributions inside each set. Literature contains methods to compare two distributions: Kullback-Leibler divergence, earth mover's distance, ... But I couldn't find anything for more than two distributions. I can imagine an indice built on top of one of these methods. For example, I could first merge all probability distributions into a single one. And then use the mean (or median) of the symmetric KL (or any other distance/divergence/...) between the merged distribution and all others distributions in the considered set. But before starting something out of the blue, I would rather be sure that nothing already exists.


title: When should one use the Reduced Major Axis regression, aka Geometric Mean Functional Relationship?

text: I recently learnt about the Reduced Major Axis regression, a.k.a Geometric Mean Functional Relationship, least products regression, diagonal regression, line of organic correlation, and least areas line (see wikipedia). The idea is that for a traditional univariate OLS regression (y=a+bx), the errors that are minimized is the "vertical" distance (along the y axis) of observations to the regression line. The reduced major axis regression assumes that there are errors to both y and x. What is minimized is the _product_ of the "y-distance" and "x-distance" of observations to the line. So it's actually the area of the triangle between the observation and the line that's minimized. Paul Samuelson apparently said it was a cool regression. However I have a hard time really understanding why and when I should use this. So, do you use it? Why?


title: How does one construct the likelihood function of a distribution in the alpha stable family given non-i.i.d. data?

text: Taking a simple alpha-stable distribution, the Normal Inverse Gaussian distribution for example, how would one derive the likelihood function _provided non-i.i.d. data_ , e.g. a price series?


title: Rare Events Model with a discrete index as the dependent variable?

text: I'm looking for a rare events model where the dependent variable is a discrete index, which means I cannot use the rare events logit model (Gary King). My dependent variable is an index of integers that range from 0 to 15, but the median is 0 and the mean is about 0.25, which implies that getting a value of 1 or more is pretty rare for this dataset. Any ideas on a model that would be better than simply running an OLS on the indexed dependent variable? And if you could recommend an R package that can get this job done, then all the better!


title: What analysis to use with a categorical DV and a within subjects design?

text: > **Possible Duplicate:**   >  Giving more details to my last post: I'm very lost and don't know which > tests to use or how to enter data in SPSS I conducted an experiment using a within subjects design and 2x2 categorical fixed factors design with a categorical DV (every subject was in every condition). I have no idea what test to use. Thanks for your help, Emily


title: Granger /Wald test interpretation

text: I have the following results from a Granger test. Can someone tell me how to interpret them? I understand that p.val(ue) < 0.05 is significant for rejecting the null hypothesis. BUT do ftest or r2(squared) add any additional information?               A g-cause B      ftest = 7.523583 p.val = 0.006165445 R2 = 0.005235041      ftest = 5.535103 p.val = 0.004030749 R2 = 0.007986544      ftest = 3.637428 p.val = 0.01242363  R2 = 0.008143064          B g-cause A      ftest = 0.3040705 p.val = 0.5814277 R2 = 0.0002147836      ftest = 0.5627182 p.val = 0.5697845 R2 = 0.0008759121      ftest = 2.365794  p.val = 0.0693209 R2 = 0.005935409


title: How to statistically test whether breastfeeding support influences breastfeeding initiation and duration?

text: I would like to know which technique is appropriate for analyzing data collected via self-administered survey. Some background information follows: Topic: Breastfeeding support and the workplace Setting: University campus Sample size: 30-40 working women who were on maternity leave from 2011-2012 (first 6 months) Hypothesis: Breastfeeding support influences breastfeeding initiation and duration. Example questions: socio and demographic items in order to describe women, Likert scale questions addressing: university policies and work culture, immediate manager support, immediate colleague support and work flow. There will be some close questions related to the intention to breastfeed. This builds on my previous question about the sampling approach and what to treat as the unit of analysis. I really appreciate all your help and consideration. Best,


title: Gaussian Process regression for high dimensional data sets

text: Just wanted to see if anyone has any experience applying Gaussian process regression (GPR) to high dimensional data sets. I'm looking into some of the various sparse GPR methods (e.g. sparse pseudo-inputs GPR) to see what could work for high dimensional data sets where ideally feature selection is part of the parameter selection process. Any suggestions on papers/code/or various methods to try is definitely appreciated. Thanks.


title: Confidence Band for Probability Mass or Density Function

text: I'm a beginner -- don't hesitate to overexplain. I need to make inferences about a population distribution from a sample distribution of a continuous variable (which is still useful as an ordinal discrete one). What I want is a confidence band on the probability function and the right method for getting the density on a given interval (or mass) with a confidence band. I thought about treating the variable as ordinal-discrete and doing proportion estimates across all the levels (taking care with the 0 "success" ones). Like I said, I'm a beginner, but I sense this is very very wrong. What's area/method should I explore? It seems like estimated density confidence bands (assuming this is what I need) are fairly recent -- I hope for something mainstream.


title: Random effect model result of R and JMP

text: I have two groups of people, namely A and B. We have the hourly mean heart rate for 45-80 hrs, length differed by individual. We are interested in the group*time effect on the heart rate. Since I have different repeats for different times in each subject, subject ID is considered a random effect. I tried lme in R and "fit model" in JMP, but get different results, the effect, such as time 10, 11, or interaction terms, such as group[A]*time[3], group[A]*time[4], are significant in JMP outputs, while not in R. Can somebody tell me if my model is correct and why are they different? In R, my script is: lme.RR<-lme(meanRR~Group*timeofday, random=~1|SubjectID) (Time of day is specified as factors) In JMP, I chose Analyze -> fit model, chose "meanRR" as Y, "group","time","group*time" in model effects, and add "subjectID" and specify it as random effects. ![enter image description here](http://i.stack.imgur.com/QzTHs.png)


title: generate frequency distribution from database

text: In a MySQL database there's a table with a single numeric `value` column. I want to plot the distribution of these values as a bar chart/histogram with the following requirements:   * there should be a maximum of N bars (intervals) in the chart   * the width of each bar (x-axis range) should be uniform and the height of each bar should reflect the number of values in this interval.   * the endpoints of the bar should occur at round numbers. I understand that this is a fairly vague requirement, but hopefully the following example will illustrate what I mean   * the intervals should be contiguous, e.g. the next interval should start where the previous one ended   * ideally it should be possible to retrieve the data with a single query   * it's OK to have intervals with a count (y-axis value) of 0 ## Example If N = 3 and the table contains the following data               +------------+     | value      |     +------------+     |     -49.2  |     |     -28.2  |     |      13.3  |     |      23.3  |     |      51.4  |     |      77.9  |     +------------+      On inspection it's fairly easy to see that the intervals `{-50..0, 0..50, 50..100}` satisfy the requirements for this dataset and value of N. However, I'm struggling to come up with a general solution that works for any value of N and any dataset. Here's what I've tried so far: ## Calculate Interval Width Get the max and min `value` via the following query               SELECT min(value), max(value), count(*) FROM my_table      Then pass the result into this (Groovy/Java) method to calculate the width of each interval               // intervalCount is what I've referred to as "N"     static Integer getRoundedIntervalSize(Double min, Double max, Integer intervalCount) {         Number intervalSize = Math.ceil((max - min) / intervalCount)              Integer roundingScale = Math.log10(intervalSize) - 1         Number roundNearest = 10 ** roundingScale              // round up the interval size to the nearest roundNearest         Number intervalDelta = roundNearest - (intervalSize % roundNearest)         intervalSize + intervalDelta     }      ## Get Frequency Distribution I then use the following query (with the value returned by `getRoundedIntervalSize` substituted for `:groupSize`) to get the number of values in each interval               SELECT      floor(value / :groupSize) * :groupSize  as groupLowerLimit,                 count(*) as groupCount     FROM        my_table     GROUP BY    groupLowerLimit     ORDER BY    groupLowerLimit ASC      This returns the lower limit of each interval and the number of values in each interval, which is all I need to build the frequency distribution. ## Shortcomings Although this approach works fairly well when the dataset is relatively uniformly distributed, when this is not the case, it results in intervals which have different widths or are not contiguous. Also, when the range of the dataset is small (e.g all values between 1 and 4) and N is large (e.g. 30), the number of intervals generated tends to be much smaller than N. Is there a better approach to solving this problem that meets the requirements above?


title: When to use regularization methods for regression?

text: In what circumstances should one consider using regularization methods (ridge, lasso or least angles regression) instead of OLS? In case this helps steer the discussion, my main interest is improving predictive accuracy.


title: Using autocorrelation to find commonly occurring signal fragments

text: I have a sensor which is capturing accelerometer data as a person walks. What I'm interested in extracting is each signal fragment when a step is taken. The Z-axis is what is used since only one axis is required to detect changes in steps. The images below illustrates a sample Z-axis gait signal (for 400 iterations).![Z-axis from accelerometer during normal gait cycle](http://i.stack.imgur.com/nZyPt.jpg) The image below illustrates the first-half of the above signal (for 200 iterations). ![Z-axis from accelerometer during normal gait cycle first half](http://i.stack.imgur.com/aCP0X.jpg) The subject is initially standing still and then begins walking at ~X=30. Notice how there is an apparent pattern as the user walks. What I'm interested in using Autocorrelation to smooth the Z-axis signal using Matlab to smooth the signal (based on the image below). Unfortunately, I don't have a strong signal processing background, and I have a decent grasp of Matlab. How can I go about achieving smoothing the gait signal so I can extract steps? The literature that I'm using suggests that steps may be extracted by looking at the peaks of the smoothed signal. ![Goal](http://i.stack.imgur.com/TaxQC.png) Other sources have suggested the use of Hidden Markov Models to extract each of the gait cycles, but I thought about a simpler signal processing approach before I consider using something advanced. However, what would be the best strategy if I wanted to pursue this course of action?


title: How can I integrate R with PHP?

text: I would like build a web application using R. I am using Windows Vista and have an Apache server. I have tried Rpad, but I was not able to correctly configure it. How do I set up Rpad as I am not that well off with PHP and Apache server? Or are there other ways to use R on Apache server?


title: How to separate distributions from weighted dataset?

text: I’m trying to separate two component distributions of an apparent finite mixture from a weighted dataset (determined by a weighted.histogram). I've a set of data and weights only for a part of the set, like this:               Data   Weights     2.5    2.8     2.6    2.8     3.4    2.8     3.2    2.8     3.5    2.8     3.7    2.8     3.1    2.8      4.4    1.0     4.5    1.0      4.5    1.0      5.2    1.0     5.4    1.0     5.4    1.0     5.3    1.0     5.5    1.0     5.6    1.0     5.7    1.0     4.7    1.0      6.2    1.0      I can separate it without the weights with the `normalmixEM` function in the `mixtools` package for `R`, but it doesn't work with weighted data. Can anybody help me to separate distributions from weighted data?


title: How to identify contributors in groups of people

text: I have a set of 120,000 people. These people are assigned to 10 million groups. Each of these groups have a target variable indicating the kind of output it made (A, B, C or D). I'd like to identify common set of people that contributed to a certain output, e.g. "groups that had Mr.X, Y and Z most likely resulted in result A". What kind of analysis is recommended to solve this case, most preferably with SPSS, R and/or Excel? I figured it would make a huge matrix of 120,000 people (columns) x 10 million groups (rows), but I wasn't sure how to move on from there. I thought of using Decision Tree, but the matrix is too big to process. (Edited on July 23) Each of these people belong to an organization (apart from the groups he/she is assigned to), so if I do a simple clustering analysis, I know I will get people clustered by the organization they belong to, and that's not the result I want. If I were to make clusters, I'd like to make clusters of people that contributed to results A, B, C and/or D. Would you have any suggestions on what kind of analysis I can perform to get such result?


title: Hyper-volume of the $\alpha$ contour of a multivariate Gaussian

text: I'm looking for the assymptotic ($n\rightarrow \infty$) value of (the log of the determinant of) the covariance of the $\alpha$% of observations with smallest Eucledian distance to the origin in a sample of size $n$ drawn from, say, a bivariate standard Gaussian. \--The hyper-volume of an ellipse is proportional to the determinant of its covariance matrix, hence the title.-- \--By standard bivariate Gaussian, I mean $\mathcal{N}_2(0_2,\pmb I_2)$ where $0_2$ is a vector of 0 of length 2 and $\pmb I_2$ is the rank 2 identity matrix.--- It's easy to see by simulations than when $\alpha=52/70$ the number is around $\approx -1.28$:               library(MASS)     n<-10000     p<-2     x<-mvrnorm(n,rep(0,p),diag(2))     h<-ceiling(0.714286*n)     p<-ncol(x)     w<-mahalanobis(x,rep(0,p),diag(p),inverted=TRUE) #These are eucledian distances, because the covariance used is the identity matrix     s<-(1:n)[order(w)][1:h]     log(det(cov(x[s,])))      but I don't recall how to obtain an exact expression (or failing that, a better approximation) for this.


title: Mean reversion in AR(1) process

text: Mean reverting level in following AR(1) process is $b/(1-a)$. $$x(t) = a + bx(t-1).$$ _I understand this._ I read that the mean reverting level for AR(1) process given below with finite differencing is $-(b/a)$: $$x(t)-x(t-1) = a+ bx(t-1).$$ Can someone explain how is this true? Can you please help me with the derivation?


title: Achieving high recall for smaller class in unbalanced linear svm

text: I have an svm-related question. I have an unbalanced dataset, meaning classA could be 1/10 to 1/35 of classB. Well I am interested in getting a _linear_ svm which would separate the data and would achieve the greatest possible recall of classA. Let me show some **examples** to be more specific:   * If 1000 samples are mapped to classA by the hyperline I want to have 80%-90% of classA. This does not mean that classA would have only 1000 samples, it could have 2000 or 10000. I am not interested in the accuracy of classA, i.e. the standard criterion to get the more samples of classA by svm out of the total samples belonging to classA.   * If I can get 100 samples of classA and no of classB it would be preferable than getting 500 samples of classA and 500 of classB. My question is how do I choose this hyperplane. I have tried the following:   1. I used unbalanced weights for the samples (inverting the class cardinalities) but I can get high classA accuracy but low recall.   2. I tries to use the default weights, i.e. for all being 1, but I got peculiar results: when run on _train set_ for all samples it gets an accuracy of 95% (as expected since dataset is unbalanced). **BUT** if checked in classA members it gives 0%!  So have you got any ideas besides playing around with weights? P.S.1 I am using liblinear (in matlab if that matters). P.S.2 I know svm optimization criterion basically does not account for what I am looking for since margin maximization maximizes accuracy but I thought I would give it a try.


title: Forgiving measure for external cluster validation

text: I'm using external validation to measure the success of a clustering algorithm. I don't consider my categories to be definite, so I'm looking for a measure that is forgiving to the following extent:   * If two clusters are merged into one, then this shouldn't be unduly penalised as it is still a good match (but still to be penalised to some extent to prevent the algorithm being pushed towards generating huge clusters)   * If one cluster is split, into two, then this shouldn't be unduly penalised   * Suppose there are two clusters, A and B. Suppose half of A is put into A2, half of B into B2 and the other half of A and B are combined into C. This kind of alternative categorisation should be penalised, more so than the first two occurrences, but not unduly as it could quite possibly represent another possible valid classification What is a good measure for this kind of cluster validation? This is related to my previous question


title: Issues related to neural network toolbox of matlab

text: I am trying to use the neural network toolbox of matlab. However, I am seeing some issues. If I use optimization methods like trainlm, trainrp, it converges well and the error is low during the training with validation set. However, when I use methods like traingd or traingdm, the error is really high. The method doesn't learn the right weights. What could be the issue?


title: penalized regression applications in epidemiology

text: I am seeking advice on penalized regression models for selecting covariates in epidemiological studies. A difficult tasks I face is feature selection while still attempting to account for confounding factors. During my academic studies, I have been taught various methods including stepwise and purposeful selection procedures; however, the later feels more like an "art" than science. I have recently learned of penalized regression models through independent studies, although applied to predictive purposes. I am wondering if there are credible sources that would suggest penalized regression for epidemiological studies to be an efficient feature selection procedure?


title: How to compare repeatability (ICC) of different groups?

text: I have calculated the ICC values for two groups and would now like to compare the ICC values to determine if the groups differ in their repeatability. In the literature people have simply used t-tests to compare repeatability but it is unclear to me how to do this. For example, with the dummy data:               ID  gr  day behaviour     1   1   1   0.361     2   1   1   0.232     3   1   1   0.240     4   1   1   0.693     5   1   1   0.483     6   1   1   0.267     7   2   1   0.180     8   2   1   0.515     9   2   1   0.485     10  2   1   0.567     11  2   1   0.000     12  2   1   0.324     1   1   2   0.055     2   1   2   0.407     3   1   2   0.422     4   1   2   0.174     5   1   2   0.613     6   1   2   0.311     7   2   2   0.631     8   2   2   0.283     9   2   2   0.512     10  2   2   0.127     11  2   2   0.000     12  2   2   0.000      I can get the repeatability measures for group 1 and 2 as follows:               library(ICC)     g1 <- ICCest(ID, behaviour, data=dummy[dummy$gr=="1",])     g2 <- ICCest(ID, behaviour, data=dummy[dummy$gr=="2",])      But how can I now determine if the repeatability of group1 is different from group2?


title: What is the probability that the i-th element is larger than all the preceding ones?

text: Homework question: code:               public int getMax(Array[]) {       int maxValue = Array[0]       int i = 1       while(i != Array.length) {         if(Array[i] > maxValue) {           maxValue = Array[i] }       }       i++     }      Question: Assume each permutation of n different integers in an array is equally probable, what is the expected or average number of times the assignment "maxValue == Array[i]" will be executed?


title: How to remove seasonality from daily electricity demand

text: I want to remove seasonality from daily electricity demand (a time series). My understanding is there is weekly (high demand on Tue, Wed, and low demand on Sat, Sun) and annual seasonality (high demand on Winter and lower on Summer). I tried to build a model to forecast daily electricity demand in R, and plot my data as shown below: ![Daily electricity demand](http://i.stack.imgur.com/pz5OS.png) I tried to remove seasonality with the following:               demand.xts.diff<-diff(demand.xts,lag=1,difference=1)     demand.xts.diff<-diff(demand.xts,lag=7,difference=1)      I also tried to use `lag=365` and `lag=366` (I am not sure what lag to use, due to the leap year issue), but none of them successfully removed seasonality. The ACF and PACF are shown below: ![ACF](http://i.stack.imgur.com/OV7T4.png) ![PACF](http://i.stack.imgur.com/HcXKi.png) Any advice is appreciated.


title: What is a standard deviation?

text: What is a standard deviation, how is it calculated and what is its use in statistics?


title: Member of multiple normally distributed populations

text: I'm writing a little tool for work and a stats question has arisen that is beyond me. In layman's terms: Let's say I've got a sample of 50 people who read the WSJ and their ages are normally distributed with a mean of 45 and a standard deviation of 10. I've got a sample of 100 people who read the NYT. Mean age: 40, std. dev. 8. And I've got a sample of 36 people who read E Weekly. Mean age: 30, std dev: 15. Then this new user comes in, and I don't know his age, but I know he reads WSJ AND NYT AND E WEEKLY. How do I quickly find the age that is most likely? Do the sample sizes matter? Thanks!


title: Where does the term "covariate" come from in statistics?

text: Wikipedia has an acceptable definition of _covariate_ , which matches exactly the Dictionary of Epidemiology. According to Webster, the first known use was in 1965. But the Google n-gram viewer has it showing up earlier, although I suspect many of the early documents are erroneously dated. None of this really answers my question, though, which is: where does this now-ubiquitous term come from?


title: Proof of Theorem 7.3 in book of Probabilistic Graphical Models by Daphne Koller

text: I'm studying graphical models myself and reading contents about bayesian networks. When I am reading in page 371, section 8.1.4 Linear-Gaussian models, in Pattern Recognition and Machine Learning, I cannot derive myself the equation of $$\mathrm{Cov}[x_i, x_j] = \sum\left(w_{jk} \mathrm{cov}[x_i,x_k]\right) + I_{ij}v_j $$ where $x_k$ are the parents of $x_j$, especially I don't know how the second term $I_{ij}v_j$ comes. I try to find the proof in other references, like in Probabilistic Graphical Models by Koller. There is a theorem 7.3 in the page 251 which states the same thing as the above equation, but it also does not supply a proof.


title: Checking the significance of association in R

text: I have 2 variables each having 500 samples. And I want to get a $p$-value which would give me the significance of association between those 2 variables. What is the way of getting this $p$-value in R? I checked `fisher.test` and `chisq.test` functions, but I think both of those require integer-valued variables while my variables are both real-valued. And when I try `cor.test` function, the p-values are quite big (about .05) even when the correlation between the variables is quite high (about .99), which does not make sense to me.


title: Bank loan model (spr?)

text: ![enter image description here](http://i.stack.imgur.com/xIFSm.png) I am trying to fit a model investigating the amount of the loan (or a transformation of it) as a function of the the variables income, gender, customer and age. Fitting a standard linear regression gives poor results (R-squared of only 0.27). Besides, the relationship between income and loan seems quite weird. It has a sort of left-rotated Z-shape (see picture) Would be appropriate to use a semiparametric model here? If so, how would you model it? I tried this in R:               model <- spm(loan~f(income, basis="trunc.poly")+age+factor(gender)+factor(customer))      But this gave poor results (no significant variables). Anyone has suggestions?


title: Finite mixture models - Basic understanding

text: I have been reading lecture slides about Dirichlet Process. ![enter image description here](http://i.stack.imgur.com/Grw1t.png) In page 22, there is a picture about the following finite mixture model. $$\phi _{k}\sim H\\\ \pi \sim Dirichlet(\alpha /K,\dots,\alpha /K)\\\ Z_{i}\rvert\pi \sim Discrete(\pi )\\\ x_{i}\rvert\phi _{z_{i}}\sim F(\cdot \vert \phi _{z_{i}})$$ **I know the meaning of the following variables (Please correct me if I were wrong):** $N:$ Number of parameters/observations $K:$ Number of mixture models $\alpha:$ Dirichlet parameter $\pi:$ Probability distribution on N variables $x_{i}$ observed variables or data that we want to model $H:$ Hyper-parameters, prior distribution about the K mixture models $\phi _{k}:$ Parameters of the k-th model $F(\cdot \vert \phi _{z_{i}})$：Individual mixture models **My questions are:**   1. What does $z_{i}$ mean and what is the relation with $\pi$ and $x_{i}$?   2. What does $Z_{i}\rvert\pi \sim Discrete(\pi )$ mean?   3. What does $x_{i}\rvert\phi _{z_{i}}\sim F(\cdot \vert \phi _{z_{i}})$ mean? Thanks!


title: Calculating joint probability problem

text: This is the example from Casella and Berger 2nd page 146. you have $$f(x,y)=e^{-y}, 0<x<y< \infty $$ and you are interested in finding $$P(X+Y\ge1)$$ They solved this problem the following way $$P(X+Y\ge1)=1-P(X+Y<1)=1-\int_0^{1/2}{\int_x^{1-x}{e^{-y}dydx}}\\\= 1-\int_0^{1/2}{(e^{-x}-e^{-(1-x)}dx)}=2e^{-1/2}-e^{-1}$$ My question is why did they solve it this way? 1) Why can't I just solve $P(X+Y\ge1)$ rather than $1-P(X+Y<1)$? I am confused because the example before this was solved different way. There $f(x,y)=6xy^2$ and to solve $P(X+Y\ge1)$, $$P(X+Y\ge1)=\int_0^1\int_{1-y}^16xy^2dxdy=9/10$$ 2) Why is the limit of integral for x is 0 to 1/2? Where do I find this value? They solve this using a graph, but is there way to justify their method without using a graph?


title: How does weka's decision table handle numerical (continuous) output and numerical inputs?

text: I've been using WEKA's decision table majority classifier in order to do a feature selection on a data set with both numerical (continuous) and categorical inputs, and numerical (continuous) output. I would like very much to know how does Weka deal with the numerical inputs and the numerical output for this case (decision table). Any other information about how the whole thing works will be very appreciated. (everything except the paper of Kohavi "The power of decision tables ..." from 1995 which I already read).


title: Standard Deviation of an Exponentially-weighted Mean

text: I wrote a simple function in Python to calculate the exponentially weighted mean:               def test():       x = [1,2,3,4,5]       alpha = 0.98       s_old = x[0]            for i in range(1, len(x)):         s = alpha * x[i] + (1- alpha) * s_old         s_old = s            return s      However, how can I calculate the corresponding SD?


title: Testing comorbidity adjusting for patient characteristics

text: I have data of patients of a large hospital with columns representing different mental and physical illnesses and patients constituting the rows. The data are in binary form (1 = patient has disease; 0 = patient does not have the disease). I want to test, given a certain mental disease, whether there are significant differences between the occurance of the physical diseases adjusted for demographic characteristics. That is, are certain physical diseases "overrepresented" in e.g. depressive patients controlling for age and gender? My supervisor wanted to calculate some kind of analysis of covariance with the mental ilnesses as the grouping factor, which seems inappropriate given the nature of the data. Which alternative tests could I use to answer the question?


title: Gaussian Process and Expectation Propagation time complexity?

text: What's the time complexity of training a Gaussian process and its Expectation Propagation approximation? (Before studying them, I'd like to understand if they are even feasible for my application)


title: Difference test for beta–binomial variables

text: Let $X \sim \textrm{Beta–Binomial}(\alpha_1, \beta_1, n)$ and $Y \sim \textrm{Beta–Binomial}(\alpha_2, \beta_2, n)$ and write $\mu_1 = \frac{\alpha_1}{\alpha_1 + \beta_1}$, $\phi_1 = \alpha_1 +\beta_1$, $\mu_2 = \frac{\alpha_2}{\alpha_2 + \beta_2}$, $\phi_2 = \alpha_2 +\beta_2$ for brevity. I have some samples $x_1, x_2, \ldots, x_k$ and $y_1, y_2, \ldots, y_k$ and want to test the null hypothesis $\mu_1 = \mu_2$. I am interested in all of the following cases, especially the third one:   1. We know $\phi = \phi_1 = \phi_2$ in advance.   2. We know that $\phi_1 = \phi_2$ but not the exact value.   3. It is possible that $\phi_1 \ne \phi_2$. I believe that in the limit $n \to \infty$, $\phi \approx n$, by approximating $X$ and $Y$ with appropriately scaled normal distributions, these situations correspond to the paired difference $Z$-test, Student's $t$-test and Welch's $t$-test, respectively. However, I am more interested in the cases when $n$ and $k$ are small, say $n = k = 10$. It would be very interesting if there was a workable solution for $n = 1$ (beta–Bernoulli distributions), but I doubt that.


title: Why (mathematically) is the parametric bootstrap usually better than the empirical one?

text: As I know from experience, the parametric bootstrap performs better in terms of coverage probability for confidence intervals then the empirical bootstrap. Of course, this makes sense because you put in some information about the distribution and that is used to reduce the error. However, thats not really a mathematical description of what is going on... Books such as Efron/Tibshirani dont go into details. Van der Vaarts Asymptotic Statistics only covers higher order consistency for studentized vs unstudentized intervals. I think Hall doesnt write anything about it either... Can i read up on this somewhere?


title: ordered logit-difference-in-difference

text: I ran an ordered logit difference-in-difference estimation and got stuck in terms of interpretation. The dependent variable takes on values 1 2 and 3 and I am trying to obtain coefficients for the variables for each of the three categories as I see in papers but I am not sure how to go about it.


title: Interpreting / understanding intercept and slope values for an interaction using summary output in R

text: In R, using lme4, I am using a mixed model to test how a response (mode; continuous and normally distributed) is affected by a covariate (species richness (sr); continuous and normally distributed) and a factor (func.group; three levels plant, herbivore, predator). I also have a random effect (community) which is just an identity. This is my model:               nr.test1<-lmer( mode ~ sr * func.group + (1|community), data=nr.test, REML="FALSE")      I plan to use model comparisons and backwards stepwise elimination of no- significant model terms, but my question regards simply how I interpret the Fixed Effects summary output below. Particularly, what is the intercept composed of? i.e. Is the intercept the y intercept (and not a slope value), and does it represent the interaction of sr with herbivores? (herbivores is the first level alphabetically) If the intercept is just func.groupherbivore, then how do I know how much different the intercept of sr:func.groupherbivore is from func.groupherbivore given that it is not listed? I also wished to know, if the intercept is indeed the y-intercept and the values below that are the differences in y-intercept for each model term, level of factor, and interaction, then how can I derive slope values to make linear (regression) equations?               Linear mixed model fit by maximum likelihood      Formula: mode ~ sr * func.group + (1 | community)         Data: nr.test         AIC   BIC logLik deviance REMLdev      62.57 68.75 -23.29    46.57   58.38     Random effects:      Groups    Name        Variance Std.Dev.      community (Intercept) 4.36963  2.09037       Residual              0.25666  0.50662      Number of obs: 16, groups: community, 6          Fixed effects:                            Estimate Std. Error t value Pr(>|t|)         (Intercept)           -0.84682    1.27221  -0.666   0.5157         sr                     0.49306    0.05798   8.504 2.44e-06 ***     func.groupplant        3.34250    0.92840   3.600   0.0047 **      func.grouppredator     6.82647    1.02830   6.639 5.08e-05 ***     sr:func.groupplant    -0.11878    0.05724  -2.075   0.0637 .       sr:func.grouppredator -0.40567    0.05019  -8.083 9.07e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Correlation of Fixed Effects:                  (Intr) sr     fnc.grppl fnc.grppr sr:fnc.grppl     sr           -0.700                                             fnc.grpplnt  -0.459  0.444                                      fnc.grpprdt  -0.482  0.473  0.754                               sr:fnc.grppl  0.293 -0.287 -0.904    -0.611                     sr:fnc.grppr  0.645 -0.831 -0.667    -0.846     0.537


title: What is the relationship between independent component analysis and factor analysis?

text: I am new to Independent Component Analysis (ICA) and have just a rudimentary understanding of the the method. It seems to me that ICA is similar to Factor Analysis (FA) with one exception: ICA assumes that the observed random variables are a linear combination of independent components/factors that are non-gaussian whereas the classical FA model assumes that the observed random variables are a linear combination of correlated, gaussian components/factors. Is the above accurate?